#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use Getopt::Long;
use Digest::MD4;
use FindBin;
use lib ("$FindBin::Bin/..");
use parse_swissprot;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<'USAGE';

USAGE:

    panda_sprot_to_fa  [--verbose] [FILE]

    Mandatory arguments:

         --dbname STRING    the postgreSQL database

    Optional arguments:

        [--verbose]         dots printed to STDERR every 2000 records
        [--user STRING]     the postgreSQL user name
        [--host STRING]     the postgreSQL database server
        [--password STRING] the postgreSQL database password

        [FILE]              The SwissProt file you wish to process e.g. swissprot.dat
                            This can be specified on the command line or via STDIN


        [--help]

RESULTS:

    STDOUT: The resulting fasta file in NCBI nr-like format with NCBI taxonomy IDs


    This script converts SWISSProt data files into FASTA format, retaining
       the taxonomical NCBI ID.

USAGE


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#    Command line options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# options variables
my $user			= 'postgres';
my $password		= '';
my $host			= '';

my $help			= '';
my $dbname			= '';
my $verbose;


#output
GetOptions(	'verbose'		=> \$verbose,
			'user=s'		=> \$user,
			'password'		=> \$password,
			'host=s'		=> \$host,
			'dbname=s'		=> \$dbname,
			'help'	   		=> \$help);

die $usage if ($help);
die "\nERROR:\n\tMissing Arguments [dbname]:\n" .$usage
				unless ($dbname);

my $md4 = new Digest::MD4;

	

	
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#	uniquify

#			removes duplicate entries in an array

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

sub uniquify(\@)
{
	my %hash;
	@hash{@{$_[0]}} = ();
	@{$_[0]} = keys %hash;	
	return @{$_[0]};
}	


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#	line_wrap

#			puts a newline every n letters
#			n.b. Do not use for constant parameters eg "$3"

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

sub line_wrap($$)
{
	my @array;
	my $seq = $_[0];
	while (length $seq)
	{
		push @array, substr $seq, 0, $_[1], "";
	}
	return join "\n", @array;
}






#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	per_record_call_back

#		parse_pfam::parse_records will call this routine once per PFAM record

#		Description
#			A) Process ID AC DE OS and OX fields corresponding to
#				the SwissProtID, gene name, description, Taxa name, and Taxa ID
#			   Print out accession and sequence in fasta format


#		Parameters:
#			1) record data: a reference to a hash whose keys are field names
#							 e.g. $pfam_data->{'DE'} is the description field
#							 these are all arrays (ie. $pfam_data->{'DE'}[0..n])
#							 except for $pfam_data->{'SQ'} which contains the sequence
#			2) file position of the record start
#			3) record number (indexed from 1) of the current record
#			4) user defined data. In this case, it is an array of
#								taxid_to_name lookups
#								file handles for pubmed, omim codes, variants

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
sub per_record_call_back($$$$)
{
	#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
	#	FASTA Sequence

	my ($pfam_data, $recstart, $reccnt, $custom_data) = @_;
	my ($taxid_to_name		,
		$obsolete_taxids	,
		$verbose) = @$custom_data;

	#print progress dots
	if ($verbose && $reccnt % 2000 == 0)
	{
		print STDERR "\n\t" if (($reccnt / 2000) % 50 == 0);
		print STDERR '.';
	}

	# make sure all required fields are present
	for (qw (DE ID OX AC SQ OS))
	{
		die "\nFatal Error:\n".
			"The field '$_' is missing ".
			"in record # $reccnt beginning on line $recstart\n"
			unless exists $pfam_data->{$_};
	}

	
	# AC
	# accession code (name)
	# only keep the first (primary) accession
	$pfam_data->{'AC'} = $pfam_data->{'AC'}[0];
	$pfam_data->{'AC'} =~ s/;.*//;
	
	
	# DE
	# join description line and remove duplicate spaces
	$pfam_data->{'DE'} = join (' ', @{$pfam_data->{'DE'}});
	$pfam_data->{'DE'} =~ s/ +/ /g;

	
	# ID
	# SP-TrEMBL ID
	die "\nFatal Error:\n".
		"More than one ID line in record #$reccnt beginning on line $recstart\n"
		unless @{$pfam_data->{'ID'}} == 1;
	# remove everything after the first whitespace
	$pfam_data->{'ID'} = $pfam_data->{'ID'}[0];
	$pfam_data->{'ID'} =~ s/[ \t]+.*//;


	# OX
	# NCBI TAX IDs
	$pfam_data->{'OX'} = join ('', @{$pfam_data->{'OX'}});
	unless ($pfam_data->{'OX'} =~ /NCBI_TaxID=\s*(\d+)/)
	{
		warn "\nWarning:\n".
			"The 'OX' line did not begin with 'NCBI_TaxID='\n".
			"in record ($pfam_data->{'AC'}) #$reccnt beginning on line $recstart\n".
			"OX   $pfam_data->{'OX'}\n";
		return;
	}
	$pfam_data->{'OX'} = [ $pfam_data->{'OX'} =~ /(\d+)/g];
	uniquify(@{$pfam_data->{'OX'}});

	# OS
	# taxonomy names, used only as backup if out of sync NCBI Tax IDs used!!
	$pfam_data->{'OS'} = '['. join ('', @{$pfam_data->{'OS'}}).']';



	# create accession line in the NCBI nr style:
	# Duplicate the whole entry for each taxid and concatenate with the ascii \x01 symbol
	my @acc_fragments;
	foreach my $tax_id (@{$pfam_data->{'OX'}})
	{
		my $acc_fragment = 	'SP|' . $pfam_data->{'AC'}.			# accession code
							'|NCBI_TAXID|'.$tax_id.'|'.			# taxid
							$pfam_data->{'ID'}.' '.				# gene name
							$pfam_data->{'DE'};					# description

		# add tax name in brackets corresponding to the taxid
		if (exists $taxid_to_name->{$tax_id})
		{
			$acc_fragment .= '['.$taxid_to_name->{$tax_id}.']'
		}
		else
		{
			$obsolete_taxids->{$tax_id}++;
			$acc_fragment .= $pfam_data->{'OS'};
		}
		push(@acc_fragments, $acc_fragment);
	}
	my $acc_line = '>'. join("\x01", @acc_fragments);

	print $acc_line,							"\n",	# accession
			line_wrap($pfam_data->{'SQ'}, 60),	"\n";	# sequence
	

}


# connect to database
my $dbh = DBI->connect( "dbi:Pg:dbname=$dbname",
                        $user,
                        $password,
						{
                          RaiseError => 0,
                          AutoCommit => 1,
						  PrintError => 0
                        }
                      )
		or die "Database connection to $dbname failed:\n$DBI::errstr\n";


		
		
		
#_________________________________________________________________________________________
#
#	Get NCBI TaxID to Taxonomy name correspondences from the Panda Database...
#
#_________________________________________________________________________________________
sub get_taxid_to_name($\%)
{
	my ($dbh, $taxid_to_taxname) = @_;
	print STDERR "\tGet NCBI TaxID to Taxonomy name correspondences from the Panda Database...\n"
		if ($verbose);
	
	((my $sql = <<"PL/SQLCMD") =~ s/^\s+/ /gm);$sql =~ s/\n/ /gm; $sql =~ s/ +/ /gm;
	SELECT tax_id, tax_name FROM taxa;
PL/SQLCMD
	my $rows = $dbh->selectall_arrayref($sql);
	$taxid_to_taxname->{$_->[0]} = $_->[1] for  (@$rows);
}





print STDERR "\tGet NCBI TaxID to Taxonomy name correspondences from the Panda Database...\n"
	if ($verbose);
my %taxid_to_taxname;
get_taxid_to_name($dbh, %taxid_to_taxname);





#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	main

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

print STDERR "\tParsing Swiss Prot file...\n\t[" if ($verbose);

my $fh_SwissProt = *STDIN;
if (@ARGV)
{
	open SP_FILE, $ARGV[0] or die "Error\n\tCould not open $ARGV[0]\n$@\n";
	$fh_SwissProt = *SP_FILE;
}

my %obsolete_taxids;
my $reccnt = parse_records(	$fh_SwissProt, 				# swiss prot file handle
							@{[qw (ID OX AC DE OS)]},	# desired sprot fields
							\&per_record_call_back,		# call back function
							[	\%taxid_to_taxname,		# has map taxid to names
								\%obsolete_taxids,		# get a list of taxids
														# without valid names
								$verbose]
								);				# ..

# swiss prot seems to use taxids which are no longer present in NCBI!
print STDERR "]\n\t$reccnt\tpfam records were parsed.\n" if ($verbose);
if ($verbose && scalar keys %obsolete_taxids)
{
	print STDERR "\t", scalar keys %obsolete_taxids, "\ttaxids were found which did not ".
				"correspond to any in the NCBI list.\n".
				"Either they are obsolete in Sprot or\n".
				"you need to rerun panda_load_ncbi_taxonomy_tree.pl with up to date files.\n";
}

print STDERR "\n";
